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Abstract 

In this paper we develop a Bayesian procedure for estimating multivariate stochastic 
volatility (MSV) using state space models. A multiplicative model based on inverted 
Wishart and multivariate singular beta distributions is proposed for the evolution of the 
volatility, and a flexible sequential volatility updating is employed. Being computationally 
fast, the resulting estimation procedure is particularly suitable for on-line forecasting. 
Three performance measures are discussed in the context of model selection: the log- 
likelihood criterion, the mean of standardized one-step forecast errors, and sequential 
Bayes factors. Finally, the proposed methods are applied to a data set comprising eight 
exchange rates vis-a-vis the US dollar. 

Some key words: multivariate time series, stochastic volatility, GARCH, state space 
models, Bayesian forecasting, Kalman filter, Wishart distribution. 

1 Introduction 

Over the last two decades, considerable effort has been devoted to the development of time- 
varying volatility models and related computational algorithms. It is widely recognized that 
volatility modeling has important implications for the analysis of returns on stocks and ex- 
change rates. More recently, attention has moved to examining the implications of volatility 
for other financial applications such as derivatives pricing, optimal portfolio selection, and risk 
management (for instance, to enable efficient forecasting of Value-at-Risk) . Although several 
univariate volatility models have been developed and are routinely used, the time-changing 
feature of the volatility is better described by multivariate models that explicitly account for 
cross-correlations among asset returns. A multivariate framework is desirable because assets 
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can be formally linked together and can be influenced by common unobserved factors; as 
a consequence of this, we often observe related movements between markets, or sectors, or 
exchange rates. 

The many efforts to model multivariate volatility fall into two main classes of models: mul- 
tivariate generalized auto-regressive heteroscedastic (M -GARCH) models an d multivariate 



stochastic volatility (MSV) models. The review paper bv lBauwens et al.l [20061 ] well describes 
the capabilities and limitations of M-GARCH models. In brief, the large number of param- 
eters, which are typically specified by maximum likelihood estimation, and the fact that the 
unobserved volatility is not modelled as a stochastic process, somehow limit the applicability 
of these models. On the other hand, MSV models are more flexible, because the volatility 
is assumed to change stochastically according to a latent process. However, most stochastic 
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Liesenfeld and Richard 



20061 ]. need essentially to resort to 



stochastic simulation schemes such as Markov chain Monte Carlo methods (MCMC), which 
may be heavily computationally intensive. Although much progress has been made on the 
front of simulation-based procedures, and more efficient algorithms are now available, the iter- 
ative nature of such procedures hampers the applicability of multivariate stochastic volatility 
estimatio n in real-time applications w here, for instance, prompt user interventions may be 



required [Salvador and Gargalld . 120041 ] . For such reasons, it would be desirable to rely on 
analytic solutions that translate into fast and flexible algorithms, while still enjoying some of 
the advantages offered by MSV models. 

Computational solutions that trade off the complexity of the model for speed are valuable, 
and have been explored in the literature. A simplification that facilitates the development of 
inferential procedures is to assume that the volatility follows a random walk (RW) evolution. 



of 


Quintana and West 




1987 


1. 


Putnam and Quintana 


199< 


H. 


Quintana and Putnam 


1996 


1, 


West and Harrison 


1997 


1, 


UhliJ 


1997 


1, 


Liu 


2000 


1. 


Sover and Tanveri 


2006 


1. 


Carvalho and West 


2007 


], and references therein. For instance 


. Harvev et al. 




1994 


] suggest an approximate in- 



ferential method for a MSV model based on the extended Kalman filter using crude mean 
and variance approximations; although the evolution of the volatility matrix is defined as an 
autoregressive (AR) process, the authors suggest that a RW evolution works equally well. 

In this work we elaborate on some of the results that have already been proposed in the lit- 
erature mentioned above. Using the co nvolution of t he Wishart and singular multivariate beta 
distributions, which was first proved inlUhlia 199411. we const r uct a RW model for the evolu 



tionof the yolatility. In the works of lAguilar and Westl 20001 ] . iLiul [2000] 



12006], and 



Carvalho and Westl [2007^ . all adopting the RW assumption, the multivariate 
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volatility estimator s resemble their counterp a rt univariate estimators b ased on gamma and 



beta distributions 



West and Harrison 



19971 . iTriantafvllopoulosl . 
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However, we have 



noticed that these estimators are incorrectly derived, in that they give rise to a shrinkage 
volatility evolution, which is not a realistic choice. In particular, we demonstrate how the 
multivariate beta density has been overlooked in the above references to the point that the 
updating equation for the degrees of freedom has been wrongly computed. The resulting 
volatility estimator proposed in this paper is a weighted average of the square logarithmic 
returns. Thus, with proper choice of the weights, the modeller obtains volatility estimators 
that guarantee mean reversion over time and are appropriate to analyze volatility. 

This paper is organized as follows. Section 12.11 defines the model and the Bayesian esti- 
mation procedure is given in Section 12. 21 Section 12.31 is concerned with model assessment and 
selection, and three performance measures are derived, namely the log-likelihood criterion, 
the mean of the standardized one-step forecast errors, and sequential Bayes factors. Section 
[3] applies our methods to a data set comprising eight foreign exchange rates vis-a-vis the US 
dollar. A proof of Section [2.31 can be found in the appendix. 



2 Stochastic volatility 
2.1 The model 

Consider a p-variate vector of log-returns {yt}t=i ... N, where t is the time index, for some 
positive integer N. The zero-drift conditional volatility model assumes 

y t = T 1 t /2 e u e t ~N p (0,I p ), t = l,...,N, (1) 

where Tt is the conditional volatility matrix of yt, €% is p-variate innovation vector following 
a p-variate Gaussian distribution with zero mean vector and identity covariance matrix; fi- 

1/2 

nally, T t den otes the square root of T t, using the Choleski decomposition or the spectral 
decomposition [Gupta and Nagarl . I200C ] . 

At time t, let y t = {yi, ■ ■ ■ ,yt} denote the information set, comprising data up to time 
t = 1, . . . , N. In order to estimate Tt, we need to define an evolution law for Tt. A sensible 
law postulates that 

E(£7 + >')=E(£ t -V), (2) 
namely the expectation from time t to t + 1 remains unchanged, and 

Var(vecp(E- + 1 1 )|y t ) > Var(vecp(S i - 1 )|y*), 

where vecp(S^~ 1 ) denotes the column stacking operator of the covariance matrix T^ 1 . These 
assumptions define a random-walk type evolution law for T^ 1 , i.e. S^Lj = T^ 1 + Tt, where 



3 



has zero mean . Suc h an evolution is possible under the multiplicative law of covariance 
matrices of 



Uhlid [19941 ] . that is 

£-\ = kUp^YBt+iUpt 1 ), t = 0, 1, 



,N-1, 



(3) 



where ^/(SJT 1 ) denotes the upper triangular matrix of the Choleski decomposition of S^ 1 , so 
that SJ -1 = Here Bt+i follows, independently of S^ -1 , the singular multivari- 

ate beta distribution (whose density is given in equation (|A-ip of the appendix). Initially, we 
assume the inverted Wishart prior 

1 



E ~IWJn + 2p,S ), 



n 



1 



(4) 



with density function 



P(S ) 



| 5o |(n +P -l)/2 etr( _ 5oS -l ; 



2P(n+p-i)/2r p ((n + p - l)/2)|S |( n+2 P)/ 2 ' 

where < 5 < 1 is a discount factor, \Sq\ is the determinant of So, etr(.) stands for the 
exponent of a trace of a matrix, and T p (.) denotes the multivariate gamma function. It is 
also assumed that the innovation sequence {e^} is uncorrelated and that {et} is uncorrelated 
with So, i.e. K(ete' s ) = (for any i / s) and E(etvecp(So)') = (for all t). From the above 
inverted Wishart prior it turns out that Sq -1 follows the Wishart distribution with n + p — 1 
degrees of freedom and scale matrix Sq 1 , i.e. Sq 1 ~ W p (n + p — 1, Sq 1 ). 

In order to completely specify this model, a value for the parameter k has to be specified. 
In Section 12.21 it is shown that in order to guarantee the expectation invariance property ([2]) 
of the RW model, it is necessary to specify k as 

5(1 —p)+p 



k 



5{2-p)+p-l 



(5) 



2.2 Estimation 



Suppose that at time t, the posterior distribution of St is 

Stly 4 ~IW p (n + 2 P ,S t ), 



(6) 



where n = 1/(1 — 5) and St is known. For the singular multivariate beta density of -Bt+i, 
we write B t+ i ~ B p (m/2,l/2), where m = 5(1 — 5)~ l + p — 1. The "singularity" of the 



distribution derives from 1 < p — 1, for any p > 1 and so the matrix I p — Bt+-\ is 



Diaz-Garcia and Gutierrez 



s mgu 



1993]). 



ar 



(for more details the reader is referred to lUhlid [19941 ] and 
The choice of m is conveniently made so that two of the assumptions of the beta density are 
satisfied, that is m > p — 1 and (1 — S)n has to be an integer (see also the last paragraph of 
Section [ 



4 



Since E t 1 \y t ~ W p (n+p — 1, S t 1 ), from the evolution ([3]) and from lUhligl 1994 ]. it follows 
that AT^E^Iy* ~ W p (n + p- l,^" 1 ) or E^Jy* ~ W p (n + p- l,fc£ t _1 ) and so the prior 
distribution of E^+i is 

E t+ i|y* ~IW p (5n + 2p,fc- 1 5 4 ). (7) 

From ([6]) we have E(E 4 ~ 1 |y t ) = (n + p — l)<5 f _1 and from ([7]) we have E^^-Jy*) = (5n + p — 
^kS^ 1 , and so by equalizing these two expectations we obtain 



A: 



n + p — 1 



5(1 - p) +p 



<5n+p-l <5(2-p)+p-l 

as in ©. Using properties of the Wishart distribution, and adopting k as proposed above, 
one can verify that Var(vecp(E^ 1 )|y*) > Var(vecp(E^~ 1 )|y t ), thus the RW type evolution ([3]) 
is verified. 

Proceeding now with the posterior distribution at time t + 1, we apply Bayes theorem 
by noting that the likelihood function from the single observation yt+i is p(yt+i|Et + i) , which 
from ([1]) is the p-variate Gaussian distribution iVp(0, Ef + i). Thus 

p(y t+ i|E t+ i, ?/)p(E m |?/) 



p(E m |y m ) 



oc 



which is proportional to 



etr(-^ +1 E- + 1 1 ^ +1 /2)|^ 1 5 f |^- 1 )/ 2 etr(-A:- 1 g t E t " + 1 1 /2) 

|E i+1 |V2|S t+1 |(^+2p)/2 
| St+1 |-(^+l+2 P )/2 etr( _ (?/t+iy / +i + k -l St )^j2), 



E i+1 |y i+1 ~/W r p (n + 2p ) 5 m ), 



(8) 



where S t+1 = fc _1 5t + y t +iy' t+1 , since 5n + 1 = n. 

Equations ([6]), ([7]) and ([8]), together with the prior (J1J) constitute a full algorithm, for 
t = 1, . . . , N— 1. We remark that, for p = 1 and /c = 1/5, the above results reduce to the usual 
algori thm for univariate stochastic volatility estimation, as reported in 



19971 ] and 



West and Harriso 



o 



Triantafvllopoulod [20071 ]. 



For p > 1, we see that, since 5 < 1, we have 5(1 — p) + p > 5(2 — p) + p — 1 and so 
< k~ l < 1. Thus by expanding as 



St = k-'So + ^ViVj > * = 1, ■ ■ ■ , # , 

3=1 



we can approximate St by 



0) 



5 



and exclude the influence of the prior Sq, which anyway is deflated as t increases. We note 
that St is just a weighted average of the log-returns {yjy'j}j=i,...,t with weights A;" 1 . Prom 
this it follows that even if {S^ -1 } follows a random walk, the estimator St is still capable 
of exploiting mean reversion of the log-returns (as it is a weighted average of the squares of 
log-returns) and thus it is a suitable estimator for the volatility. The posterior mean of £j 
and the prior mean at £t+i can be derived easily from the inverted Wishart densities, i.e. 



E(Et|y* 



S t 



[I " S)S t 



and E(S m |y*) 



k^St (1 - 5)S t 



n-2 25-1 " v t_ri| " ' 5n - 2 k(35 - 2) 
the posterior mean being defined for 5 > 1/2 a nd the prior mean being defined for 5 > 2/3. 



In related work, a number of authors such as 



1997 . Chapter 16]. lAguilar and Westl [20001 ] 



Quintana and West 



Liu 



200C ], 



19871, 



West and 



Sover and Tanveri 



larrison 



20061 ] . and 



Carvalho and West 



20071 ] have suggested to use k = 1/5. Although it is easily verified that this is a correct choice 



when p = 1, setting k = 1/6 when p > 1 results in a shrinkage- type evolution for £ t . This 
can be seen by first noting that, with k = 1/5, we have 

E(S- 1 1 | ?/ t ) - E(£7/ V) = (p - lXcT 1 - 1)S^ (10) 

and therefore the expectation is not preserved from time t to t + 1, as we have E^^ly*) > 

In particular, when p is large, even if 5 ~ 1, the above model postulates that the estimate of 
£^l_i is larger than that of Sj -1 - In other words {EJ -1 } follows an AR model EJT 1 = aE^ + r^, 
where a > 1; such a setting is clearly inappropriate. With the RW type evolution of ST/ 1 , 
claimed in all the above references, assuming that the limit of St exists, it follows from (|10p 
that = (p — l)(5~ l — 1) lim^oo St. This, for p > 1, implies that 5 = 1 or lim^oo S^ 1 = 0, 
two meaningless results. Our suggestion is that 5 should be replaced by as in (JSJ), a 
choice that now preserves the expectations. 

Furthermore, for p > 1, the updating equation of the degrees of freedom of the Wishart 
distribution suggested in the above references, namely 

nt + 2p = Snt-i + l + 2p = n <5* + (1 - <5*)/(l - 5) + 2p, 

does not seem to be correct. The reason for this lies in the multivariate singul ar beta distri- 
bution, Bp(m\/2 ) m2/2) which is only defined for be ing a positive integer Uhlid. I1994H . 



West and Harrison 



19971 ] and ISover and Tanveril [20061 ] , re- 



Setting mi = (1 — 5)nt, as in 
suits in m2 not being a positive integer. In our algorithm, we resolve this issue by setting 
nt = n = 1/(1 — 5) so that mi = (1 — 5) n = 1. For more details on the multivariate sin gular 
beta distribution the reader is referred to 



Srivastava 
appendix. 



Uhlie 



|l994 



Dfaz-Garcfa and Gutierre2 



19971 ]. and 



20031 ]; the density function of this distribution is given in equation (|A-ip of the 



6 



2.3 Performance measures 



2.3.1 The likelihood function 

One method of model judgement and model comparison is via the likelihood function. In this 
section, first we derive the likelihood of our model in closed form. Adopting approximation 
0, the only parameters that need to be selected in order to fully specify the model is the 
scalar 5, since k is specified in ([5]). Using the following result of Theorem [IJ one possibility 
is to choose the value of 5 that maximizes the log-likelihood function (under the restriction 
2/3 < 5 < 1). 



Theorem 1. In model ([2P-(Ep the log-likelihood function o/Ei, . . . , Ejv, based on data y%, 



is 



N x N N 



35-2 
2(i~5) 



N 



E lo sl s *l' 



t=l 



for 



Np N Np(25 - 1 

c = log 7r log 2ir ; — 

2 B 2 6 2(1 -5) 



io g * + ivio g . r P {2-Hi-srHHi-P) + P)} 



r p {2-m-s)-i(s(2-y) + y-i)y 

where 5 > 2/3, k is as in |5|j and Lt is the diagonal matrix with diagonal elements the positive 
eigenvalues of I p - k- 1 {U(E^} 1 )'}- 1 T£ 1 {U(T^ 1 )}- 1 , with S^ 1 = W(E^ 1 )'W(E^ 1 ). 

The proof of this result can be found in the appendix. A common modelling strategy 
in Bayesian inference is to plug the posterior mean of in to the likeli hood function and 
then to compare models by comparing their likelihood functions (e.g. see Leonard and Hsu] 



1999] ). Th i s app r oach has common roots to estimation methods using the profile likelihood 



Liitkepohll . 12003 . 



Leonard and Hsu 



19991 ] . and clearly it has the advantage of combining 



Bayes estimation with likelihood-based inference. In addition to that, this approach can 
be very useful for choosing nuisance parameters, such as the discount factor 5. The max- 
imization of the log-likelihood function with respect to 5 may be slow because this is a 
non-linear function in 5. A possibility would be to evaluate the log-likelihood function only 
on a few admissible values for 5 (2/3 < 5 < 1). Values of 5 lower than 0.7 can result in 
very volatile, not smooth, and thus unstable posterior estimates of E*; values of 5 larger 
than 0.95 can result in very smooth estimates of not able to capture the clusters and 
the spikes of the volatility. In this paper (see the illustration of Section [3]) , we recomm end 



exploring values of 5 in the range .7, 0.75, 0.8, 0.85, 0.9, 0.95. IWest and Harrisonl [19971 ] and 
Triantafvllopoulos and Nasonl 20071 ] have some discussion on the performance of the posterior 
estimates at the boundary values of discount factors 5 > 0.95. 
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2.3.2 One-step forecast error 

Other than the log-likelihood function, the mean of square standardized one-step forecast 
error vector (MSSE) provides another performance measure. From ([I]) the one-step forecast 
distribution of yt+i \y t is a p-variate Student t density with 5/(1 — 5) degree s of freedom, mean 
vecto r and scale matrix fc" 1 ^, written y t+ i|y* ~ t p (5/(l — 5), 0, k^St) [Gupta and Nagaii . 
2000 ]. It then follows that, for 5 > 2/3, 



Var(y t+ i|y* 



fc- 1 ^ (i - S)S t 



5/(1 -5) -2 (35-2)*;' 
which also can be derived from Section 12.21 using conditional expectations, i.e. 

Vax(y t+1 \y t ) = E(Var(y t+1 |E t+1 , y^y 1 ) = E(E t+1 \y t ) = L'^ffi , 

since from model ([1]), it is Var(E(yt+i\^t+i,y t )\y t ) — 0. Having obtained an expression for 
the variance, we can now write the standardized one-step forecast error vector ut+i as 

6 



_ x = y/kS t 1/2 yt+i with u t+ i|y* ~ t p (j-^i °> !p 



(11) 



so that the vector 



(1 - S)S t 



-1/2 



yt+i 



m !(35-2)fe 

has E(u* +1 |y') = and E(n* +1 (tij +1 )'|y*) = I p . Then the MSSE vector is given by 



1 N 

MSSE = -]T{(^) 2 ,...,(^) 2 }\ 



t=i 



where u\ - 

(i,...,iy. 



It) ■ ■ ■ 5 U pt) ■ 



Models that fit well the data are expected to yield MSSE 



2.3.3 Bayes factors 

A third approach for model d i agnos t ics is based on sequ e ntial Bayes factors 



1997 



Salvador and Gargalkl . I20041 . iTriantafyllopoulod . 120061 ], Suppose we have two compet- 



West and Harrison 



ing models, M\ and M 2 , parameterized in terms of 5\ and 5 2 , respectively. First, a Bayes 
factor is obtained as the logarithm of the ratio between the density of 114 = ut(5\) (under M\) 
and the density of u t = Ut(52) (under M 2 ). Specifically, at each time t we have 



H t = log 



t = l,...,N, 



p(ut(h)\y t ~ 1 ,M 2 ) 
and, from the Student t density (jlip . this becomes 

r((m + P )/2)r(n 2 /2) r (1 + u t (5 2 y Ut (5 2 )r^ \ 1/2 



H 



T((n 2 +p)/2)r(ni/2) I (1 + u t (5i)'u t (5i))^+P 
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Figure 1: Daily observations on eight foreign exchange rates. 



where T(.) denotes the gamma function and n« = <5j/(l — <5j), for i = 1,2. 

A value of Ht > then suggests that model M\ has to be preferred over M2, in the sense 
that A4\ is associated with a superior forecast distribution. Alternative, negative values for 
Ht suggest that M2 is the preferred model. In situation where Ht = 0, both models are 
deemed equivalent. One point of interest is what decision can we make when Ht fluctuates 
around zero. In such a case one may sele ct a threshold value in order to decide which model 
to choose, as in lWest and Harrison! 1997 ]. 



3 An illustration using foreign exchange rates 

In this section we present an analysis of eight exchange rates vis-a-vis the US dollar. The 
exchange rates are the Australian dollar (AUS), British pounds (GBP), Canadian dollar 
(CAD), German Deutschmark (GDM), Dutch guilder (DUG), French frank (FRF), Japanese 
yen (JPY) and Swiss franc (SWF), all expressed as number of units of the foreign currency 
per US dollar. The sample period runs from 2 January 1980 until 31 December 1997, and 
corresponds to 4774 observations, sampled at daily frequencies. This data set was originally 
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Figure 2: Sequential Bayes factor Ht of the standardized one-step forecast errors of model 
Mi (5 = 0.7) vs model M 2 {5 = 0.95). 



obtain ed from the New York Federal Reserve, and then discussed in iFranses and van Dijk 



20001 ] . Figure [T] illustrates the daily observations on the level of all eight exchange rates. 



We have applied the stochastic volatility model of Section [2.21 to the logarithmic returns, 



which have been collected i n a vector y± = (y 



it, 



ies of exchange rates, as in iQuintana and Westl 19871 ] . iPutnam and Quintanal 19941 ] . and 



Quintana and Putnam 



UR tY . Following the empirical stud- 



19961 ] . we adopt the random walk for the evolution of the volatility 
and thus we specify A; as in (|5]). In order to choose a suitable value for the parameter 5, we 
have used the performance measures described in Section 12.31 Following suggestions in that 



section, we have only considered a few selected values of 6 in the range 0.7 < 5 < 0.95. The 
results from this analysis are summarized in Table HJ which provides the mean of the MSSE 
(MMSSE), the log-likelihood function (evaluated at the posterior mean of the volatility), and 
the mean of the Bayes factors of the standardized one-step forecast errors. For the computa- 
tion of the Bayes factors here, each Mi is based on the current value of 5, and is compared 
against a baseline model M.2 that uses 5 = 0.95. 

From Table [Tl it can be observed that for small values of 5, the MMSSE also attains small 
values, indicating poor performance, when compared to an ideal MMSSE value of one. This 
result seems to suggest that the forecast covariance matrix of yt has been over-estimated. As 
5 gets close to one, the MMSSE also gets close to one, which underlines an improvement in 
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Table 1: Mean (over the eight exchange rates) of the mean square one-step forecast standard- 
ized errors (MMSSE), log-likelihood function (LogL) evaluated at the posterior mean of the 
volatility, and mean of the log Bayes factor H t (t = 1, . . . , 4774). 



5 


MMSSE 


LogL 


H 


0.70 


0.072 


-12857.59 


-6.269 


0.75 


0.194 


-12395.30 


-5.681 


0.80 


0.337 


-11721.93 


-4.950 


0.85 


0.506 


-10644.32 


-3.982 


0.90 


0.701 


-8627.03 


-2.564 


0.95 


0.912 


-3458.23 






the estimation of the forecast covariance matrix of yt- The log-likelihood function attains its 
largest value at 5 = 0.95. For each 5 < 0.95, the Bayes factor mean H is negative and this 
indicates a preference in favour of model M.2 (for 5 = 0.95). In particular we note that the 
model performance deteriorates as 5 decreases, a fact that is captured by all three diagnostic 
measures considered here. As a result of this, we conclude that 5 = 0.95 produces the best 
model. 

Figure [2] shows the log-Bayes factor sequence {Ht}, from which the superiority of model 
M.2 is clear. We observe that, out of N = 4774 data points, {Ht} is positive at only 37 
points (i.e. only 0.77% of the time). Using sequential Bayes factors, the modeler has the 
extra advantage of choosing the discount factor at each time t according to the sign of Ht- 
This is particularly advantageous in an on-line setting, and when decisions have to be made 
in real time. 

Figure O shows the posterior volatilities, i.e. the estimates of ou : t (i = 1,...,8), for a 
subset of the data points (t = 4001, ...,4774). Most of the volatilities are small, except for the 
JPY/USD; even for small volatilities, this figure indicates clearly the highly volatile periods 
for each exchange rate. Figure H] shows the posterior correlations of GBP/USD versus all the 
other rates. This figure confirms that the correlations are time-varying. By inspecting Figure 
S]we observe that GBP/USD is most correlated with DUG/USD, FRF/USD, GDM/USD, 
and SWF/USD. 

Finally we note that, for this relatively large data set, based on 4774 time points in 8 
dimensions, the estimation algorithm (implemented in the R language on a Windows plat- 
form) took less than a minute (55 seconds) to complete, on a PC with Intel(R) Celeron(R)M 
Processor 1.60GHz and 504MB RAM, including the evaluation of the log- likelihood function 
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Figure 3: Estimate of the posterior volatility for the FX data, using the model with 5 = 0.95. 



and the Bayes factors. 



4 Conclusions 



In this paper we have described a Bayesian modeling approach for multivariate stochastic 
volatility. The proposed estimation methodology is delivered in closed form, is easily imple- 
mentable and efficient, as the model relies on only one parameter. 

The models proposed in thi s pape r are closely related to the above m entioned articles as 
well as to the models of 



Uhligi 19971 ] and iPhilipov and Glickman 



20061 ] . Notably, we have 



shown that similar volatility estimators proposed in the literature are based on a shrinkage- 
type volatility evolution, which is not a realistic choice. Instead, the estimator described here 
guarantees a random walk type evolution. 

The procedure proposed in this paper attempts to combine the simplicity of non-iterative 
algorithms with the sophistication of stochastic volatility models. In our view, algorithms such 
as the one suggested here are particularly attractive because they can model high dimensional 
data with low computational cost, which is crucial for certain real-time applications in modern 



12 




if 
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1000 2000 3000 4000 
(b) GBP/USD and CAD/USD 




1000 2000 3000 4000 
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1000 2000 3000 4000 
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1000 2000 3000 4000 
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1000 2000 3000 4000 
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1000 2000 3000 4000 
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Figure 4: Estimate of the posterior correlation coefficient for the exchange rates data using 
the model with 5 = 0.95. 

computational finance, such as algorithmic trading. Future research efforts will be directed 
towards other financial applications with special focus on optimal portfolio allocation. 



Appendix 

Proof of Theorem^ First we derive the density of E^Ef—i, for t 
B t ~ B p (m/2, 1/2), for m = 5(1 - +p - 1, with density 



1, . . . , N. From ©, it is 



p(B t 



7T 



-p/2 



r p ((m + l)/2) 



if f |-p/ 2 |a|( m -P- 1 )/ 2 



(A-l) 



where /„ 



r p (m/2) 

Bt = HiKfH[, Kt is the diagonal matrix with diagonal elements the positive 



eigenvalues of I p — Bt, and H\ is a matrix with orthogonal columns, i.e. HiH[ = I p . For 



more details on this distribution see 



Uhlie 



199 



Now from evolutio n fl3l) we have the transformation from Bt to St = k 1 (W(S t _ 1 1 )) 1 B t 1 



x(W(S t l 1 1 )')" 1 . From 
is 



Diaz-Garcia and Gutierre2 



19971 ] the Jacobian of this transformation 



dB t ) = |Kt| p / 2 |Ltr p/2 AT p/2 |St-i| 1/2 (dS 



t), 
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where Lt is the diagonal matrix with diagonal elements the positive eigenvalues of I p — 
/c~ 1 (^(S t l 1 1 )') _1 S i T 1 (^(S f l 1 1 )) _1 . From the above transformation it is 

|W(S-\)| = \U(E-} 1 )'U(E-\)\ 1 / 2 = \Z t -i\ 1/2 and \B t \ = ^S^HEtl" 1 
and thus from (|A-1[) 

= ^ r ^^+ i y 2 ) fc -^| S ,_ 1 |V2 

1 p\frij 

x\K t \ ~p/ 2 | S t | (— P-i)/ 2 1 tf t |P/2 1 Lt I -p/ 2 

= 7r"P/ 2 /t"P( m -p)/ 2 lKi!!i±llZ^l|L l-p/2 

r p (m/2) 1 " 

x |S t _i | (™-p)/ 2 | S t |-(^-p-i)/2 _ ( A _ 2 ) 

For the likelihood function L(E; y), where E = (Si, . . . , Sat) and y = {y\, . . . , yN), write 

v 

L(S;y) = JJp(yt|E t )p(E t |S t _i). 

t=i 

From equation ([I]) we have yt|St ~ N p (0, St), while the density of Sf|Sf_i is given by (|A-2j) . 
The required formula of the log-likelihood function is obtained by taking the logarithm of 
L(S;y). □ 
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